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CHAPTER 1 


INTRODUCTION 

The scattering of elastic waves by an inhomogeneity embedded in an 
infinite homogeneous isotropic elastic medium has been studied by 
numerous investigators. An "inhomogeneity" is a region in which dif- 
ferent material properties from its surrounding medium exist. Eshelby 
[1,2] developed the method of equivalent inclusion to determine the 
elastic field of an ellipsoidal inclusion. An "inclusion" is considered 
to be a region which has the same geometric shape and dimension as the 
inhomogeneity but the same material properties as its surrounding medium 
after an eigenstrain is imposed within that region. Mai and Knopoff 
[3] appeared to be the first in applying Eshelby' s result to form the 
scattering theory of a single sphere. Subsequently, Gubernatis [4] 
also used Eshelby 's result to study the long-wave scattering of elastic 
waves for an ellipsoidal inhomogeneity. In the above studies Eshelby' s 
solution for the static displacements was used as the first approxima- 
tion in the iteration. 

Wheeler and Mura [S] first applied the concept of eigenstrain to 
study composite materials, in which they considered the difference in 
elastic moduli between the inhomogeneity and matrix. Subsequently, 

Fu [6] presented a formulation for the elastodynamics field of two 
ellipsoidal inhomogeneities embedded in an infinite elastic medium sub- 
jected to plane time-harmonic waves. He [7] later gave a complete 
formulation in extending the method of equivalent inclusion to dynamic 
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elasticity and gave some results for three-layered and five-layered 
media subjected to plane time -harmonic longitudinal waves. The scat- 
tering of plane waves by an ellipsoidal inhomogeneity is presented in 
[8,9]. 

The scattering of a plane compressional wave by a spherical inhomo- 
geneity in an infinite elastic medium has been studied by Ying and 
Truell [10] and by Pao and Mow [11]. They used the method of separation 
of variables to solve the wave equation which describes the incident, 
reflected (scattered) and refracted waves inside and outside the spheri- 
cal inhomogeneity. Because the solutions are expressed by a spherical 
coordinate and the inside and outside scattering field can match exactly 
along the boundary of the spherical inhomogeneity, the results by this 
method are considered as exact solutions. Some numerical results of 
the elastic scattering cross section were shown by Johnson and Truell 
[12], and some dynamical stress concentration factors around a spheri- 
cal cavity were found by Pao and Mow [14] . 

This research is concerned with numerical calculations according 
to the extended method of equivalent inclusion by Fu [6-9] , and a com- 
parison with the exact solutions. Two cases are studied here. The 
first case is concerned with the three -layered problem and the other 
case is concerned with the spherical inhomogeneity problem. Results 
by the direct integration method for the three-layered problem and the 
method of separation of variables for the spherical inhomogeneity 
problem are compared with those obtained by the extended method of 
equivalent inclusion applied on the two cases. 
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CHAPTER 2 

WAVE MOTION IN AN ELASTIC MEDIUM 
2.1. Governing Equation 

Consider the infinitesimal element of an elastic body with mass 
density p in the absence of body forces, the equations of motion are 

*jk,k " P 5 j 

where a dot indicates a differentiation with respect to time while a 
subscript comma indicates spatial differentiation. 

If the elastic material is linear and homogeneous, the stress- 
strain relation is 

a., = C., e (2) 

jk jkrs rs 

where C.. are the elastic constant, 
jkrs 

If the elastic material is isotropic, then the independent elastic 
constants are reduced to two and 


C jkrs ’ ls 3 k 5 rs * »( s jr s ks * 5 js 

where A, p are Lame's constants and <5^ are kronecker's delta. 

Substituting eqs. (2 . 1 . 2) (2 . 1 . 3) into eqs. (2.1.1) the equations 
of motion in terms of the displacements are obtained as follows: 


(A+y) U. .. + u U. .. = p U. (4) 

J,Ji i,JJ i 


or, in vector notation as: 


(5) 


(A+u) VV • if + = p f 
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where V is the vector differential operator. 

The displacement vector can be decomposed as 


tf=Vt(; + ^x<£, V • tp = 0 . (6) 

where the ip and $ are the scalar and vector displacement potentials, 
respectively. 

A substitution of eqs. (2.1.6) into eqs. (2.1.5) leads to 


V[(A+2ti)V 2 i|) - p tjj] + V x[y V 2 ^ - p £] = 0 


(7) 


and the 41 and 41 satisfy the wave equations 



( 8 ) 

(9) 


2 2 

in which v^ = (A+ 2 y)/p and v^ = y/p are the velocities of longitudinal 
waves and shear waves, respectively. 


2.2. Initial and Boundary Conditions 

In the above section the general wave equations in the absence of 
body forces have been mentioned without considering external forces. 
Usually the external force can be treated as an internal or surface 
source (incident wave) , or in the form of initial conditions (impulses) 
and boundary conditions (displacement or force or mixed boundary condi- 
tions) . Here the external force from sources, i.e. incident waves, will 
be discussed. 

If the incident waves are assumed to be plane waves, in general 
they are composed of compressional and shear waves. When the incident 



5 


waves move in an elastic medium, the congressional and shear waves will 
propagate independently. But they can not travel independently if there 
is an inhomogeneity in the elastic medium. When the incident waves im- 
pinge on the inhomogeneity, compressional and shear waves will be re- 
flected back into the matrix while the same types of waves will be re- 
fracted into the inhomogeneity. Both the reflected and refracted waves 
must satisfy the general wave equations. If the elastic inhomogeneity 
is bounded to the matrix at all times, then the tractions and displace- 
ments must be continuous at the interface between the inhomogeneity and 
matrix . 


2.3. Time-harmonic Wave Problem 

If the incident waves are time -harmonic, then the reflected and re- 
fracted waves will also be time-harmonic waves of the same angular fre- 
quency. Therefore, the displacement of all the waves can be represented 
by 

U(r,t) = U(r) exp(-iwt)e u (1) 

where to is the angular frequency and e u is the unit vector in the direc- 
tion of wave propagation. 

The scalar and vector displacement potentials are 


*(r,t) = <H?) 

exp (-itot) 

(2) 

$(?,t) = <f>(?) 

exp (-icot)e^ 

(3) 


where e is the unit vector in the direction of shear wave. 

4 > 

Substituting eqs. (2 . 3 . 2) (2 . 3 . 3) into eqs . (2.1.8)(2.1.9), the 
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wave equations are reduced to 


(V 2 + a 2 ) 

\|j(t) =0 

(4) 

(V 2 + 8 2 ) 

<f>(?) =0 

(5) 


where a=w/v^, 8=aj/v^, are the wavenumber of longitudinal and shear waves. 


respectively. 
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CHAPTER 3 

INHOMOGENEITY PROBLEMS 

3.1. Transmission and reflection of waves in an infinite three-layered 
medium 

A plane compressional incident wave which is simple harmonic is 
assumed to propagate in a three-layered medium. The geometry and ma- 
terial properties of the medium are shown in Fig. 1. The displacement 
of the incident wave is 

0^ = U Q exp(ik_*r_ - io>t)e z (1) 

where U Q is the amplitude and k, u are wavenumber and angular frequency, 
respectively. e z is the unit vector in the positive direction of z-axis, 
i.e. plane wave propagation. 

3.1.1. Integration of the governing differential equation 

When the plane compressional incident wave impinges on the inter- 
face between the elastic inhomogeneity (medium II) and its surrounding 
medium (medium I), a compressional wave is reflected back into the 
medium I, while a compressional wave is transmitted into region III [13]. 
For convenience of the ensuing discussion, the displacements and stresses 
associated with the incident, reflected, refracted and transmitted waves 
will be designated by the superscript (i), (r), (f) and (t) . Details 
for incident compressional plane waves along the +z axis are given below. 
The displacement of these waves are 


Z ^ = U Q exp(ia 1 z - iwt) 

(1) 

frl 

z = -Aexp(-ia 1 z - iwt) 

/ — 
CM 

V / 


I 


8 


(to 

z 

= B exp(ia^z - iwt) 

(3) 

(f) 

z 

= U^fz) exp(-ia)t) 

(4) 


where U q , A, B are the amplitude of the incident, reflected and trans- 

2 2 2 2 

mitted waves, respectively, = ai /v^ = /(A^+2 y^), where the sub- 

script 1 denotes the matrix. 

ff y 

The displacement function U fz) should satisfy the wave equation 

(V 2 + a 2 2 )U(z) = 0 (5) 

2 

2 22 P 2 W 

where a-, = to /v T = = — t- and the subscript 2 refers to the inhomo- 

2 L U 2 +2y 2 J 

geneity . 

ff ■) 

For the simple case chosen, U (z) are obtained as 

U^(z) = Ccosa 2 z + Dsina 2 z (6) 

where C, D are constants. 

The stresses associated with the above displacements are 


zz 

= ia 1 (A 1 +2]i 1 ) U q exp (ict^z - iojt) 

(?) 

a Cr) 

ZZ 

= -iaj (A 1 +2y 1 ) A exp(-i a;L z - icot) 

(8) 

a (t) 

ZZ 

= ia 1 (A 1 +2.y 1 ) B expUc^z - icut) 

(9) 

ZZ 

= a 2 (A 2 +2y 2 3 (-Csina 2 z + Dcosct^z) 

(10) 

The stresses 

and displacements must be continuous at the 

inter- 

Thus at z = 

- 6/2 the continuity conditions are 



u^+ u Cr) = 

(11) 


z z z 


(i) + 00 = (f) 

a zz CT zz CT zz 

(12) 
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where f = (X2 + 2 y 2 ^/ and h = P 2 /Pj ’ 

The displacements and stresses in the whole medium can be obtained 
by using equations (1-4) and equations (7-10), respectively. 

3.1.2. Extended Method of Equivalent Inclusion 

The equations of motion for the inhomogeneity and the matrix are: 


jk,k = P 1 u j 

matrix 

(1) 

jk,k = p 2 U j 

inhomogeneity 

(2) 


I 
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The stress-strain relations are: 


„ _ r cu c 

jk “ jkrs £ rs 


matrix 


a = C?jP s 

jk jkrs rs 


inhomogeneity 


(3) 

(4) 


where pp and p 2 , are mass density and the elastic constants 

for the matrix and the inhomogeneity, respectively. 

Replacing the inhomogeneity with an equivalent inclusion which has 
the same material properties as the surrounding matrix after an eigen- 
strain is imposed in the inclusion, the governing equations are equations 
(1,3) and 


°jk,k ■ *1 u j 

inclusion 

(5) 

a = 

jk jkrs e rs 

inclusion 

(6) 

e * 


(7) 

e„„ = e - e 

rs rs rs 



e * 

where £ rg , £ rs , £ rs are the total strain, elastic strain and eigen- 
strain, respectively. 

The equivalence conditions derived and given in [7,9] are: 


AC jkrs * C jkrs S3 1 ’ * ' AC jkrs inside the inclusi ° n 

Ap + C?^ s = - Ap inside the inclusion 

j jkrs rs,k j 


where AC = 

jkrs jkrs jkrs 


Ap 


C™ 

jkrs 


£* v = ^ 

rs,k jkrs rs,k jkrs 


U. 

J 



= p 2 " P 1 


(II) 


'rs,k 


+ U. 
3 


(m) 


and 


( 8 ) 

(9) 


( 10 ) 

(ID 
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where the superscripts "i" and "a" denote field that is induced by the 
presence of mis-match and that is applied. 

Equations (3.1,2.5-3.1.2.11) are the general equations for the 
equivalent inclusion with no restrictions on the geometric shape of the 
inclusion, i.e. these equations can be applied to an inclusion of arbi- 
trary geometry. 

For the three-layered problem, the equivalence conditions are re- 
duced to 

Apu 2 l C* n [0]A n + ^ 1} [0]B n ] + (Xj_ + 2v 1 )A q = -Apa> 2 H o 
n=o 

TT APU) 2 l [^ 1) [0]A n + ^ 2) [0]B ] + (Xj + 2P 1 )A 1 = -Apu\ (12) 

^ " n=o 




n=o 




-Apw 2 H, 


m 


and 


jr(AX + 2Ap) l [<^ 1} [0]A n+ ^ 2) [0]B n ] + (A 1+ 2y 1 )B o = -(AA + 2Ay)] 


n=o 


li-(AX + 2Ap) l [^ 2) V^n 33 [0]B n ] + (V 2tl l )B l = -(^ + 2Ap)E 1 


n=o 


i-j-(AA+2 Ay) “ [0] A n+ <^ m+2) [0] B n ] + (X^Zy^ = -(AX + 2Ay)E m 


n=o 


( 13 ) 


where: ♦„(*) - f s ( 2 ' )V“1 1 ^ Uz- 


£ 2 
*00,^ _ d* 

^ (Z) = —jj. * (2) 


AX = X 2 - Xi 


Ap = p 2 - y 1 


(iO m 

H = i — U 

m m! o 


E = 
m 


CiajD 


m+1 


m: 


The eqs . (3.1.2.12) and (3.1.2.13) give 2 (m+1) simultaneous equa- 
tions, having matrix dimensions of 2 (m+1) x 2 (m+1) if matrix represen- 
tation is used. 

The displacement and strain fields inside and outside the equiva- 
lent inclusion are: 


U (z,t) = U (z) exp(-iojt) 


m 


■ l * VA' Z >] * V 101 ’ 

n=o 


£ zz Cz,t) = e 77! (z) exp(-iu)t) 


zz 


m 


£ zz Cz) = l [A n 4>;(z) + B n ^(z)3 + ia 1 U o e ia l Z 


n=o 


and the stress fields are: 


a zz (z,t) = ct Z2 Cz) exp(-iwt) 


a zz (z) = (Xj + 2v : ) e zz (z) 


matrix 


o,,(z) = (X, + 2 p 2 ) e zz Cz) 


zz 


(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 


inclusion 
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Evaluation and integration of the $ - integrals are shown in 
Appendix I. 

3.2. Scattering of waves by a spherical inhomogeneity 

A plane compress ional incident wave which is simple harmonic is 
assumed to move in the positive z-axis in an infinite elastic medium 
where an elastic spherical inhomogeneity is embedded. The configuration 
of the whole domain is shown in Fig. 2. The displacement of the inci- 
dent wave is 

= U Q exp(iaz - iuit)e z Cl) 

where U Q is the amplitude and a, w are the compressional wavenumber 
and angular frequency, respectively, and e^ is the unit vector in the 
positive direction of the z-axis. 

In this section two methods, the separation of variables ap- 
proaches and the extended equivalent inclusion method, are used to 
determine the displacements and stresses of the reflected (scattered) 
waves. Then, two measurable physical quantities, differential scat- 
tering cross section and the total scattering cross section far from 
the inhomogeneity, will be expressed in terms of the scattered asymp- 
totic values. The differential scattering cross section dP(o))/dQ [4] 
is a measure of the fraction of incident power scattered into a parti- 
cular direction, where dfi is the differential element of solid angle. 

The total scattering cross section P (u)) is the ratio of the average 
power flux scattered into all directions to the average intensity of 


the incident fields. 



3.2.1. Method of separation of variables 


In this section, the work by Ying and Truell [10] and that by 
Pao and Mow [11] will be briefly introduced. 

When a plane compressional incident wave impinges on the surface 
of the elastic inhomogeneity, scattering occurs. Both compressional 
and shear waves are reflected back into the matrix (Medium I) while 
the same types of waves are refracted into the inhomogeneity (Medium II) . 
The potentials, displacements and stresses associated with the incident, 
reflected and refracted waves are denoted by the superscript (i) , (r) 
and (f) . The wave equations in terms of potentials are 


0 ?“ 


(V“ 


where 


«1 2 ) 

\p ^ 

= 0 



= 0 

ai 2 ) 


= 0 

O 

^ r) 

= 0 

a 2 2) 


= 0 

8 2 2) 


= 0 

’ 8 1 

2 

= p : 

L t u 2 /y 


ctj = y = P^to /^1 3 a 2 = P 2 a) ^ (^ 2 + ^^ 2 ^ * 


(i) 


( 2 ) 


(3) 


e 2 2, 

B 2 = p 2 oj /v 2 


The displacement of the incident wave is 

= U Q exp(ia^z - iojt) 


(4) 
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and the associated potential functions are 

U 

^ UJ = ^ exp(ia 1 z) 

= 0 

can further be represented in terms of a spherical coordinate 
function 

f 1 U 

= ~ l (2n+l) i n j n Ca ir ) P n (cos0) 

1 n=o 


(5) 

( 6 ) 

(7) 


where j n (x) are the spherical Bessel functions of the first kind, and 
P n (x) are Legendre polynomials. 

Because of axisymmetry of this problem, displacements, stresses 
and potential functions are independent of the spherical coordinate <5 . 
Therefore, after solving the wave equations, the potentials of the re- 
flected wave are [12] : 

00 

= l C-i) n+1 a(2n + l)A n h^r^Ccose) (8) 

n=o 


$ (r:) = l (-i) n+1 a(2n+l)B n h (B r)P (cos0) (9] 

n=o 

where h n (x) are the spherical Hankel functions of the first kind. 

For the refracted wave inside the inhomogeneity the potentials are 

[ 12 ]: 

iK f;) = l (-i) n+1 a(2n+l)C n j n (a 2 r)P n (cos0) (10) 


4> 


Cf) = 


i t-D "* 1 


n=o 


( 11 ) 


a(2n+l)D n j n (S 2 r)P n (cos0) 
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The resultant waves in the medium I and medium II are 


<P : = ip 


= 4> 


iP 2 = P 


^2 " + 


U = |i + i(D a <()) 
r 3r r v 9 VJ 


n = k ii JL/r 
U e r 36 36 ^r^ 


2 3iJj 


rr 


ct 90 = 2p [- 


2 r r 3r 2 9" 

r 

3 2 - 2 ct 2 , 1 3 ^ 1 

2 v r 3r r , 

c°te 3 rn l n 3$, 

~r 36 (D r^ “ 7 °e 77 ] 


(i) 

+ * cr) 



(12) 

(i) 

+ 



(13) 

(f) 




(14) 

(f) 




(15) 

of 

displacements 

and stresses 

in terms 

can 

be found in 

[15]: 






(16) 





(17) 

1 

+ — 
r 

d q (~ 

Q l dr 

- £) 
r 


(18) 

V 

cot9 

2 

r 

dp 

Jd 


(19) 


a a = u4^[~ - -^=- ib + $ <p +-•£*- + + — D q i1-] 

r9 39 L r 3r 2 y y r 3r 2/ 2 9* J 

r “ 


,2, . 2 : 3$ ^ 2_ , . 2 
r r 
D q 4> 


0 r 8^-2a 2 1 dip cot9 3iJ; “9 T cot9 3 

= 2UC_ "1 * + 737 + — 70 + — - — 39 °r^ 


r 3r 2 36 

r 


( 20 ) 


( 21 ) 


19 13 3 

where D r = ^ -g^r and D g = jjne" 30’^ ■ For the incident field, 

<Jj = p ^ , <j)=0, a=Oj, 8=8^. For the reflected field, p = p^^ , 6 = <p( r ) , 


a=a^j 8=82 • For the refracted field, p=p 


Cf) .-.(f) 


, a=a 2 , 8=S 2 • 
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Consider the boundary conditions at r=a 


yCi) + jjCr) 

r r 


u' 15 * U e tr:i 


rr rr 


^ + a C ^ 


. U C£) 


= a 


= a 


r9 


r9 


r 

Cf) 

0 

Cf) 

'rr 

r0 


( 22 ) 

(23) 

(24) 

(25) 


Substituting eqs. (7-11) and eqs . (16-21) into the above continuity 
equations will form four simultaneous equations. After solving these 
equations, the unknowns A n , B n , C , D n will be obtained. Thus, the 
elastic field inside and outside the inhomogeneity can be determined. 

The total scattering cross section far from the inhomogeneity is 

[ 12 ]: 

P(“) = 4it&2 l (2n+l) [ [ A | 2 + n (n+1) -i- |B n | 2 ] (26) 

n=o P 1 n 

and P (to) can be normalized as 

00 a 

p(“) = 4 l (2n+l) [ | A | 2 + n(n+l) -i | B | 2 ] (27) 

n=o n Sj n 


3.2.2. Extended Method of Equivalent Inclusion 

The eqs. (3. 1.2.1 - 3.1.2.11) are general equations by this method 
and can be applied to the inhomogeneity problem of arbitrary geometry 
[7-9]. The equivalence conditions are recorded here as follows: 


k Apaj2 [fs-j^A-j + f s - jk [°]A ;jk +...+F cW [<)]B v , + F^.jOjB^ , + ...] 


si 


skj 


T r skjA LUJD kjr 


+A S * - Apw H s 


( 1 ) 
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rrAP» 2 [£ 5j _ p tO]A. * f sjk , p ro]V...*F skjip [0]B kj * F skjl>p [0]B kjV ...] 


+A = -Apw H 
sp sp 


( 2 ) 


* 3T 2A "[ d stj tOUj * d stjk tO]A. k *. . •*D stjk [0]B jk *D st . ki [0]B. k .*. . .] 
+ <Vst B ™* 2u l B st> * -C“ 5 st E mm * 24 >* B st> C3) 


1! A ^ 5 st ^ d mj ,p ^ A j +d mmjk,p ^ A jk + ‘ ‘ ' +D mmjk,p ^ B jk +D mmij k,p ^ B i jk + ' ' 
*fr 24 “t d st j ,pt°] A j* d 3t j k,p[ 0 ] A jk*---* D stjk, p [ () ] B jk* i:i stjki,pt») B j ki + -'-! 


+ (At 5 B + 2y..B ) = - (AA6 .E + 2AyE „ ) 

1 st mmp 1 stp J ^ st mmp stp 


C4) 


where: Ap = - p^ , AA = A^ - A^ , Ay = U 2 - 

H = U s=3 

s o 

= 0 otherwise 

H = icOJ s=3 p=3 

sp o r 

= 0 otherwise 

E = iaU m=3 n=3 

mn 0 

= 0 otherwise 

E = -U a 2 m=3 n=3 p=3 

ranp 0 r 


0 


otherwise 
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4 irp l 0 ) 2 £ js [?] = - [S l^ 5 js + *,js " 

4 irp l“ 2 f jsk [?] = -t 6 l 2 Vjs + *k,js " *k,js ] 

[r] = -tW^kj + 2 p l & l\k 6 sj 
- 2 V,skj + 2 y l*,kjs ] 

4irp l“ 2F skj£^ = -^l a l 2 ^,s 5 kj + 2y i e i 2<!> ii,k 6 sj 
- 2p lVskj + 2p l^,kjs^ 


4.7T0 ai^d frl = - + d> 5 . + 4* • — 4* • I 

4 p l" mni L J °1 L ,n jm *,m jn ,jmn ,jmn J 


4irp l w d mnjk [?] = ' S 1 C *k,n 5 jm + ^.m^n + \,jmn ' 


jmn^ 


2 

4irp 1 aj D mn T k rrl = 2u A , . l-p, 8 2 j> . 6, +<)> . 6, ] 

1 mnj k (.rj ^ kjmn- v , kjnrn J M 1 1 L ,jn km Y ,mj kn J 

- Vl^.mn^k 

4irp l w D mnijk^ = 2p l ^k, ijmn'^k, ijmn^ " P 1 6 1 ^k, in' 5 jm +<f> k,mi <S jn^ 

- X l a l\,nm 6 ij 


where the <j>- and ^-integrals and their derivatives are evaluated for a 

sphere by using the method given in [8] and are listed in Appendix II. 

From Appendix II, f s;jk [0], F s j k [°]- f sj,p [01> F skj £,p ’ d stj [0] ’ 

D s tjki [°]. d stjk,p [0] md D stjk,p t01 are equal t0 ZSr °- ^ 6qS - 
(3. 2. 2.1 - 3. 2. 2. 4) can therefore be much simplified. 

If one assumes a two term expansion in the eigenstrains, then 

eqs . (3. 2. 2.1) and eqs. (3. 2. 2.4) will form 21 simultaneous equations 

while eqs. (3. 2. 2. 2) and eqs. (3. 2. 2. 3) will form another 15 simultaneous 
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equations. The two sets of simultaneous equations are uncoupled. 

Solving the two set simultaneous equations for A^ , A^ k , B^ k and 
B , the displacements of the scattered wave can be written as [8,9] 

" £Vj + Wjk + WV + ' 

where: 4ir Pl « 2 l aj = - + # >mj - * >mj 

4irp l“ 2 ^Bji = " 8 l 2 *l 6 mj + V«j ’ 

^l^mkj = ’t X l a l 2 ^m 5 kj + 2y l 0 l 2 ^k 6 mj 
- 2 V,mkj + 2 V,mkjl 

4irp l^mkjil = + 2y l B l 2 *£,k 5 mj 

' 2w l**,®kj + 2w l*A,mkjl 


where : 


hr] = III 


exp(ia 7 R] 

t ^Vi 


r outside £2 


exp (ia 1 R) 

= II! a *i — dv ' 


r outside £2 


expCia.R) ^ 

♦let...,® - III R dv ’ r 0UtSide 




’,kl. ..s^ " 3x ]( 3x 1 . ..3x s 


♦ [?] ■ III 


exp(i8,R) 
± A\n 


r outside £2 


exp (ig.R) ^ 

Kih = HI x£ dV’ r outside £2 
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Y k& . . . s 


[?] = Iff x^x^.-.x 


expfiBj^R) 


n 


dV' 


r outside & 


> _-V 3 . r"^-i 

■*,*« = *w 


w - 


',k£ . . . s L " J 3x k 3x £ ... x s 


$[r] 


As | r ] -*- « the far field value of U m can be expressed as [9]: 

expia 1 r expiS.r 

U = C — + D 

m m r m r 


( 6 ) 


and the associated far field stress are 

expic^r 

a = iX,a, C, l, 

mn 11 r k k mn 


expia^r 
‘1 r 


iu. Fa, — (C 2. + C l ) 

*1 L 1 r m n n nr 


(7) 


expi6 1 r 

+ 0. — (D l + D & )] 

1 r m n n m J 


where are direction cosines. 

The differential scattering cross section is therefore [9] 

i 2 

ill n i i i i 

dP(aQ 
dSl 


c 

2 a, 

D 

in 

1 

m 

u 

0. 

U 

0 

1 

0 


( 8 ) 


and the total scattering cross section is 


P(ui) = / 


dn 


(9) 


and P(oj) can be normalized as 


’(“>) » ~T / 


ira 


dP(ai) 

dG 


dG 


( 10 ) 


Some important volume integral calculations are shown in Appendix II. 
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CHAPTER 4 

COMPARISON OF COMPUTATIONAL RESULTS 


In this chapter numerical- results are presented and compared for 
the three-layered and spherical inhomogeneity problems. Numerical 
results were made known by Truell and his co-workers for a perfect 
spherical inhomogeneity. In order to compare with these results, the 
same material properties are used and listed in Table 1. 


Table 1. Material Properties 


Material 

Compressional 
wave velocity 
(m/s) 

shear wave 
velocity 
(m/s) 

Mass density 
(g/ cm3) 

Stainless Steel 

5790 

3100 

7.90 

Mg 

5770 

3050 

1.74 

A1 

6568 

3149 

2.70 

Ge 

5285 

3376 

5.36 

Polyethylene 

1950 

540 

0.90 

Be 

12890 

8880 

1.87 


When the extended equivalent inclusion method is used, there are 
complex matrices, with dimensions depending upon the dimensionless 
wavenumber and the differences in elastic constant and mass density 
between the inhomogeneity and the matrix, need to be solved. In some 
cases, the dimension of the matrix are so large that the numerical 
error can be relatively high. In order to get the best results, the 
IMSL (International Mathematical and Statistical Libraries) subroutine 
LEQ2C is used to solve matrices with double precision. The routine 
applies iterative improvement until the solution is accurate to machine 
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precision. If the matrix is too ill-conditioned to get effective itera- 
tive improvement, a terminal error is produced. 

4.1. Three -Layered Medium Problem 

The input parameters for both methods are the dimensionless wave- 
number ct^6, the relative ratio of elastic constants f and mass density 

h, where f is (X 2 +2u 2 )/ (Xj+ 2^) and h is p 2 /p^. T ^ ie displacement ampli- 
tude U q and the stress amplitude (X^+2u-^)ct^U o given in the preceding 
figures are nondimensionalized. 

The calculation of the exact solutions is simple and the dimension 
of the matrix is just 4x4. The calculation of the equivalent inclusion 
method is relatively complicated. The solution is an infinite series 
summation, and N^the accepted number of terms to get convergent results) 
depends on o^S, f and h. In the calculation of the displacements and 
stresses from Eqs. (3.1.2.14 - 3.1.2.20), the summation is considered 
to be acceptable until the ratio of the current term to the current 
partial term is less than 0.5%. 

In order to make detailed comparison, three cases are studied 
here. The first case considers the difference in elastic constants, 

i. e. f^l and h=l. Fig. 3-7 and Fig. 8-12 display the displacement 
amplitude and stress amplitude vs a^5, respectively. The second case 
considers the difference in mass density, i.e. f=l and h^l. Fig. 13-17 
and Fig. 18-22 display the displacement amplitude and stress amplitude 
vs c^S. The last case considers the difference both in elastic 
constants and mass density. Fig. 23-27 and Fig. 28-32 display the 
displacement amplitude and stress amplitude vs a- 5. From the results, 
it is found that the extended method of equivalent inclusion gives 



excellent results which can be treated as the exact solutions. The 
values of in getting the convergent displacement and stress ampli- 
tude in the third case are listed in Table 2. 


Table 2. The Value of N£ for the Three-Layered Problem 


a l 6 

Ge in A1 
f =1.285 
h=l . 985 

A1 in Ge 
f=0 . 778 
h=0. 504 

Polyethylene in Be 
f=90. 790 
h=2 . 078 

Mg in St 
f=0 . 219 
h=0 . 220 

St in Mg 
f=4 . 572 
h=4. 540 

0.5 

6 

4 

6 

6 

6 

1.0 

6 

6 

6 

6 

8 

2.0 

8 

6 

6 

8 

8 

4.0 

12 

8 

8 

10 

12 

6.0 * 

16 

10 

8 

14 

14 

8.0 

18 

14 

8 

16 

18 

10.0 

22 

16 

10 

20 

20 


4.2. Spherical Inhomogeneity Problem 

The input parameters are the dimensionless wave number u^a, the com- 
pressional and shear wave velocity and the mass density of the inhomo- 
geneity and the matrix. 

For the method of separation of variables, the results of the scat- 
tering cross section were plotted in Truell's paper but the specific 
values for different a^a are not listed."*" 

For the equivalent inclusion method, there are two independent 
series summation. The first one is to get the inside function value 
in the eqs. (3. 2. 2.1 - 3. 2. 2. 4), which are related to a^a and 
only. The second one is to get the outside function value in the eqs. 

(3. 2. 2. 5 - 3.2.2.10), which also depends on a^a and 0^/a^ . 

From the three-layered problem, it is known that the value of 

tThe data shown in this section are calculated from the computer program 
kindly supplied by Dr. J. Gubematis. 
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depends on ^6, f and h. This is also true for the current problem. 
Besides, it is found that the value of also depends on It 

should be noted that the dimension of the matrix is proportional to the 
value of N£. Therefore, for large N^, the dimension of the matrix will 
be very large and the derivation and numerical computation to get the 
scattering cross section is very lengthy and time-consuming. Instead 
of finding the scattering cross section for large ct^a, it is hoped to 
find the accepted value in what range of a^a when 11^=1 . Later, the 
result is called one-term solution when is equal to 1. The eqs. 

(3. 2. 2.1 - 3. 2. 2.4) are reduced to: 


A,» 2 [f 3:j [0]A. * F skj [0]B kj ] ♦ A s ■ -A»» 2 H s 

4XS s t['W°] A :i * * 2AM[d stj [0]A j * D stjk [0]B jk ] 


+ (X-6 . B + 2y.B J = -(A AS ^E + 2AyE J 
M st mm 1 st- 7 k st mm st' 


CD 


( 2 ) 


From the formulas in Appendix II, it is found that F^j [0] and 
d . [0] are equal to zero, which make the eqs. (4. 2. 2.1) and (4. 2 . 2 . 2 ) 

Su J 

uncoupled, i.e. the difference in mass density and the difference in 
elastic constants will have no coupled effects. After some manipula- 
tion, it is found that the nonzero variables are A^, B^, B 22 anc * ®33 
with B^^ equal to B 22 b 7 symmetry. For low c^a, the closed form 
solutions of these variables can be obtained and the closed form of the 
scattering cross section can also be obtained. 

Fig. 33-36 display the scattering cross section vs a^a from two 
method. It is found if ct^a is less than 1, the tendency of the one 
term solution is good compared to the exact solution. For certain 



material system the one-term solution represents a good approximation 
up to medium frequency range. The accuracy, however, is not good 
enough to produce relative errors less than 1% though a^a goes down to 
0.01. The Fig. 37-40 display the comparison for the a^a between 0.01 
and 0.1. It is therefore necessary to take more terms in the series 
for higher ct^a if this approach is to be used. 
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CHAPTER 5 

DISCUSSION AND CONCLUSION 

It is worth studying the numerical characteristics of the ex- 
tended method of equivalent inclusion. By doing this, we will know 
what kind of inclusion this method can be applied to. Because of the 
simplicity and low computer time, the three-layered problem is chosen 
first for studying. Fig. 41-50 display the displacement and stress 
amplitude vs c^6, f and h, respectively. In these figures, the curves 
both for N (number of terms) = 6,12 and for the exact solution are 
shown. From these figures, it is found N^ will increase while a^5 
or h increases or f decreases. Besides, the value of N^ for convergent 
stress amplitude is larger than the value for convergent displacement 
amplitude. It is also found that the displacement amplitude is less 
than the exact solution and the stress amplitude is larger than the 
exact solution at the larger a^5 if N is less than N_^. 

Though the acceptable numerical results of the equivalent inclusion 
method in spherical inhomogeneity are not obtained in this study, the 
one-term solution supplies a very good tendency at low wavenumber. It 
is expected when the two-term solution (N=2) are applied, the accuracy 
at low dimensionless wavenumber will be very good and it will also 
supply good tendency in higher dimensionless wavenumber. And if N^ 
is large enough for different inclusion, the value from this method 
should be very close to the exact solution as in the three-layered 
problem. 
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From the above study, it is found will be large as the a^6 
(or a^a) or h or is very large or f is very small. Therefore, 

the dimension of the complex matrices will also be very large, which 
will probably cause the numerical singularity (algorthmic singularity) 
in high a^a range. The numerical singularity may be overcome after 
some numerical improvement skills are applied or other numerical methods 
are used. It is suggested that numerical scheme be carefully studied 
for large h, S^/ct^ or small f at high wavenumber range. 
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EVALUATION OF SOME INTEGRALS FOR THE THREE-LAYERED PROBLEM 

1.1. f h Cz') n e iot l z - z ’ldz- 

-h 

, ,,n iaz„ r , 

= (-1) e R n (ah) f or z ■> h 

. [ C -1)V“ * e" ic,z ]e lcJl R ii (oh) - I[1 - (-D r ] nlZ "' r - 

r= o (n-r) ! (ia] r+i 

for -h _< z _< h 

-iaz_ , . . 

6 RjjCaJO for z _< -h 

»h„r e : R (ch) - e iah ? (-nW- r 

IT “ 1 

r=o (n-r) ! (ia) r 1 


■ , n . , n-r 

- (-l) n e~ lah l n! h 


r=o (n-r)! (ia) r+ ^ 


n+2 5/2 


1.2. a n+ \(z) -|U-/ Cz'] n e ia i Z - z ’Uz' 


-5/2 


ia 6 


_ .. -2 r C-l) n e iaz +e~ iaz ~ . , .. 
— i [ z - 6 Lj} (ct5 ) 


- 2 - 2 k -2 


£ l+(-l) m n! (az) n 3 

m=o 2 (n-m) ! (ictz) m 


P 

where: 1 ^( 016 ) = 1 


■ a5, n 
n! ^ 


n C^-) 


m=o (n-m) ! (- 


T , n-k +2 (k ) f . n-k +2 d K 

1.3. a <j> (z) = a — r- 4 (z) 

n dz K n 


-T<- z <-T 


_ .k -2 r (-l) n e laZ +(-l) k e” iaz T*. , 

= 1 L- — z =r — e Ln(ao) 

n-k < f • > m 1 f \ n 

T lllzH nl (<*z) 1 

“ 9 n+k ^ 

m=o (n-m-k) ! (iaz) 



30 1.4. a 1 " 1 " 2 * 00 [0] 


z=0 


(1) both n and k axe not odd integer 


.n-k-2^00 [0] , i 


ictd 
.k-2 2 


e L n (ouS) 


icd 


. k — 2 r 2 , . . n ! ^ 

= x [e L n (crf ) - — ] 

i 


(2) both n and k are odd integer 

icsS 

an " k+ \ Ck) [0] = " ik " 2 ® 2 L n Ca6) 


icd 

= -i k ' 2 [e 2 L n Ca6) + % 

i 


(3) only one of n and k is odd integer 

a n " k+2 <C n ( ' k ' ) [°] = 0 


_ . -.n laz 

T c n+2 A frl - C -1 ) g- 

1.5. a ^(zJ _ 


2i 


72 - R n^^ 


_ e 


-laz 


2i 


2 n 


R„ CcsS ) 


„ , ,,n...k laz 

n-k+2 _ (-1) (i) e .. 

1.6. a <}> n (z) - 2 04 


2i 

f-n k e' iaz 

1 1J ■% R (cd ) 

2i 2 


icsS ,-cd >n 

~2~ n n ! 1 -T -J 
where R (cd } = e I - 

m=o (n-m) ! ( — J 


icd /•£?L'> n 

-t-x,-.-*- \ ^ 2 > 


, . , ,icd .m 

m=o (n-m) ! (— j ~ ) 


for n < k 


for n > k 


for n < k 


for n > k 


, 5 

for z >_ j 

for z <_ - 


for z > 


for z < - 


K)|Oj NJl 0- * 
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II. 1. 


where 


II. 2. 


EVALUATION OF SOME VOLUME INTEGRALS FOR 
THE SPHERICAL INHOMOGENEITY PROBLEM 


From Ref. [6,8], we obtain: 

«£> = /// JiEiiaSI dv , 

r< R 

n n- a 

- Ill 


t inside ft 


(-1] l k n-il-k 


■ L n 0 L n Jl!k!(n-S,-k) 
n=o l=o k=o ' 


x y z 


fff 


0 3x ,£ 3y' k 3 z ' ( - n ~*' _k - ) 


gx P i ” r ' dx* dy- dz' 


• l 


,n 


n=o 


(-1) 

- j— XX ... X 

n! p q v 


m=o 


nth power of x 
m 2m , , 

J a • y i j 'J- q '.“J -i 

(2m)! pq...v 1 ^ (2m-I) ! pq. . .v J 


00 r 2m , , 2m- 1 , . 

r V (- 1 ) a „(n),m . 7 C-l) a g (n),m 

<■ L (7m 1 I l nn_ v 1 L 


R = r - r r 


f *\ . n / i \ 2 in “ 1 

C n , m = fff p(r') — - — ^ dV' 

pq...v 3X-3X-...3X; 

s (n),m _ fff 3 n (r>) 2m - 2 dv , 

pq...v JJJ Pi - T J 3x'3x’...3x' 

r 0 p q v 

P(r') = 1 


^a...s' 


(?) = fff x^ X* ex £ l iaR I dV’ 


a 


? (-.i) n 

= y - — — xx . . .x 
L n! p q v 
n=o ^ n 


nth power of x 
.m 2m 


r inside 0 


m-1 2m-l 


r r v Z 1 i a p , • r J a c^ ik J 3 111 1 

1 (2m)! pq...v (2m-l) ! pq...v‘ 


, (n) ,m 


m=o 



where : 


r (n) ,m 
pq. . .v 


cCn) ,m 
pq. . .v 


-/// PCX') 
ft 

= /// PCx*) 
ft 


3 n (r') 2m ~ 1 

3x^3x^..3x; 


3 n (r') 2m ~ 2 

§~x’3x'...3x' 

p q v 


dV» 


dV' 


PC?’) = x’x’...x; 


In what follows we develop these formula specially for a sphere of 
radius a: 
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pq. . .v M ^ J 
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= 1 and ft is sphere 


2m+2 
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a 

2m+l 
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pq 
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pq 
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^-(2m-l)a 2m 5 

3 pq 
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•=— (2m-2)a 5 

3 pq 
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2m- 2 4n 
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15 of indices 


0 


otherwise 
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f4} _ 2(1)^111 _ 2^33,111 __ q 

P Pqu p pqu 


II. 4. C Cn) ' m and S^’ 111 for p(r') = x! and fl is sphere 
” pq . . . v 1 r 


(13 C 


pqu. . .v 
(1) ,m 


P 

(1), m 


(2m-l) 2m+2 4n „ 

2m+2 a 3 6 pi 


2m- 2 2m+l 411 
2m+l a 3 pi 


C23 c (3:),m = 

pqu 


0 for m=0 

(2m-13 (2m-33a 2m — ■ p=q=u=i 


(2m-l) (2m-3)a 2m y| 


for any two equal pairs of 
indices 


otherwise 


s (3) ,m = C 2m-2) (2m-4)a 
pqu 


2m- 1 


4n 

5 


_ P=q=u=i 


= (2m-23 (2m-43 a 


2m- 1 4n for any two equal pairs of 

15 indices 


0 


otherwise 


(3) C 


(53 ,m 
pquvw 


0 for m=0 , 1 

(2m-13 (2m-33 (2m-532m a 2m " 2 • D 


D = 

4n 

7 

i=P= 

: q=u=v=w 

D = 

4n 

r i=p 

1 

q=u=v=w 

35 

ip=q 

i=u=v=w 

D = 

4n 

105 

i=P 

q=u v=w 

D = 

0 

otherwise 


p.q^^.w can change the 
order arbitrarily 



34 


s ( 5 ),m = (-2 m _2) (2m-4) C2m-6) C2m-l)a 2m " 3 • D 
pquvw 


(4) C 


Co) ,m 


D - 4n 


7 

i=p=q=u=v=w 

D = M 

■ i=p q=u=v=w 

35 

[ p=q i=u=v=w 

D - 4n 


105 

i=p q=u v=w 

O 

II 

a 

otherwise 

S (°) ,rn _ 

c C 2 ),m = s (2) 


pq pq 


pquv 


pquv 


= 0 


II. 5 


« , 2m r v 50 , .^m-l 2m -1 , . 

a) ♦[«]- Z ^kn-cc°>.”,i Z ( - 13 f2 J, 

m=o K J m=l 


C 2m- 1 ) ! 


2m 00 f . V m — I 2m — 1 

C2> *.„t0] - (-D 2 [ Z ( (^)! Cp^' * 1 Z ‘ ,l.„l S ™ 5 ’ 1 


pq 1 


m=o 


m=l 


(2m-l) ! 


pq 


m 2m 


, m-1 2m- 1 


pquv 


m=o 


(4),m 

(2m) ! '"pquv m f ^ (2m-l) ! '"pquv 


a — r t ii* xiu r » f > iu— x 

C3) ^ n ,J0] - C-1D 4 C I 4 &f - c l ; ffl + 1 l - W- 1 


] 


C4) *,p[0] = *, pqu [0] = * Jpquvw [0] = 0 


II. 6. 


k?) 


fff e gL(.i.?.Rl dv * 
a 


R 

n n-£ 


'Ill 

n=o £=o k=o 


(-D 


£ ! k ! (n-£-k) 


3n( 12E_-£) 

2 £~ k~ n-£-k 
d x 3y 3 z 


r outside R 


j7/(x')\ y') k (z') n_£ ‘ k dx'dy'dz ' 


for | r | -> °° and R is sphere 
00 n n-£ 

|(r) = l l l F (n, £,k, a, a, r) 

n=o £=o,2 k=o,2 



where: F(n,£,k,a,ct,r) = (-ia) 


35 


II. 7. 
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-EXACT SOLUTION 

o CURRENT SOLUTION AT Z/6=-0.5 


g a CURRENT SOLUTION AT Z/ 5 = 0.5 



Fig. 3 Displacement amplitude as a function of 04 <5 for the 
three-layered problem, Ge in Al, with h=1.0 
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-EXRCT SOLUTION 

©CURRENT SOLUTION RT Z/<$ =-0.5 
CURRENT SOLUTION RT Z/<5 = 0.5 


o 
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Fig. T Displacement amplitude as a function of ^ 5 for the 
three-layered problem. Stainless Steel in Mg, 
vith h=1.0 
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-EXACT SOLUTION 

©CURRENT SOLUTION AT 1 / 5 =-0.5 
o CURRENT SOLUTION AT 1 / <5 = 0.5 



Fig, 11 Stress amplitude as a function of <5 for the three- 
layered problem. Mg in Stainless Steel, vith h=1.0 























-EXRCT SOLUTION 

©CURRENT SOLUTION RT 1 / « =-0.5 


§ a CURRENT SOLUTION RT 1 / « = 0.5 



Fig. 17 Displacement amplitude as a function of a 16 for the 
three-layered problem. Stainless Steel in Mg, 
vith f=l. 0 
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-EXACT SOLUTION 

oCURRENT SOLUTION AT Z/«*-0.5 


O a CURRENT SOLUTION AT Z/ 5 = 0.5 

O 



Fig. 21 Stress amplitude as a function of aj <$ for the three- 
layered problem. Mg in Stainless Steel, vith f*1.0 





- EXACT SOLUTION 

©CURRENT SOLUTION RT Z/5 =-0.5 


O aCURRENT SOLUTION RT Z/ 6 = 0.5 
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Fig. 23 Displacement amplitude as a function of c*i <5 for the 
three-layered problem, Ge in A1 
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-EXACT SOLUTION 

oCURRENT SOLUTION AT 2 / =-0.5 


O a CURRENT SOLUTION AT Z/ 6 =0.5 
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Fig. 27 Displacement amplitude as a function ofaj 5 
three-layered problem. Stainless Steel in Mg 


for the 
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-EXACT SOLUTION 

oCURRENT SOLUTION AT Z/ a =-0.5 


o aCURRENT SOLUTION AT 1/6 = 0.5 



Fig. 29 Stress amplitude as a function of q 5 for the 
three-layered problem, A1 in Ge 
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Fig. 3 1 * Scattering cross section as a function of <-^a 

for the spherical inhomogeneity problem, A1 in Ge 
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-EXACT SOLUTION 
oCURRENT SOLUTION N =6 


O a CURRENT SOLUTION N =12 
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Stress amplitude at the point z/£=0.5 as a function of f 
for the three-layered problem at the wavenumber alrf =2.0 
and h=l . 0 
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-EXACT SOLUTION 
©CURRENT’ SOLUTION N = S' Vl ' 


O a CURRENT SOLUTION N =12 
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Fig. 50 Stress amplitude at the point z/«T=»0.5 as a function 
of h for the three-layered problem at the vavenumber 
o|,£»10.0 and f=1.0 
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